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Abstract 

In the last decade, the study of cosmic microwave background (CMB) data has become one of 
the most powerful tools to study and understand the Universe. More precisely, measuring the 
CMB power spectrum leads to the estimation of most cosmological parameters. Nevertheless, 
accessing such precious physical information requires extracting several different astrophysical 
components from the data. Recovering those astrophysical sources (CMB, Sunyaev-Zel'dovich 
clusters, galactic dust) thus amounts to a component separation problem which has already led to 
an intense activity in the field of CMB studies. In this paper, we introduce a new sparsity-based 
component separation method coined Generalized Morphological Component Analysis (GMCA) . 
The GMCA approach is formulated in a Bayesian maximum a posteriori (MAP) framework. 
Numerical results show that this new source recovery technique performs well compared to 
state-of-the-art component separation methods already applied to CMB data. 

Key words: Blind component separation, Sparse overcomplete representations, Sparsity, Cosmic 
microwave background, Sunyaev-Zel'dovich, Morphological component analysis. Morphological diversity 
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Introduction Investigating Cosmic Microwave Background (CMB) data is of huge sci- 
entific importance as it improves our knowledge of the Universe [8] . Indeed, most cosmo- 
logical parameters can be derived from the study of CMB data. In the last decade several 
experiments (Archeops, Boomerang, Maxima, WMAP - [1]) have already provided large 
amounts of data and astrophysical information. The forthcoming Planck ESA mission will 
provide new accurate data requiring effective data analysis tools. More precisely, recover- 
ing useful scientific information requires disentangling in the CMB data the contribution 
of several astrophysical components namely CMB itself, Galactic emissions from dust 
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and synchrotron, Sunyaev-Zel'dovich (SZ) clusters [12] to name a few. In the frequency 
range used for CMB observations [3], the observed data combines contributions from 
distinct astrophysical components the recovery of which falls in the frame of component 
separation. 

Following a standard practice in the field of component or source separation, which has 
physical grounds here, the observed sky is modeled as a linear mixture of statistically 
independent components. The observation with detector i is then a noisy linear mixture 
of n independent sources {sj}j=i,-- ,n '■ Xi = Yl^=i ^ij^j + The coefficient reflects 
the emission law of source sj in the frequency band of the i-th sensor; Ui models instru- 
mental noise. When m sensors provide observations at different frequencies, this linear 
mixture model can be rewritten in a more convenient matrix formulation : 



X = AS + N (1) 



where X is the m x t data matrix the rows of which are the observed data maps in each 
channel, A is the m x n mixing matrix, S is the n x t source matrix the rows of which 
are the sources Sj, and N is the mxt noise matrix. In practice, both the sources S and 
their emission laws A may be unknown or only partly known. A component separation 
technique then aims at estimating both S and A from the data X. This problem refers 
to Blind Source Separation (ESS). 

Amongst all the physical components mixed in the observed data, each one raises sci- 
entific interest. Thus it would be worthwhile to devise a separation technique able to 
differentiate effectively between most physical components. Up to now, several source 
separation techniques have already been used in the field of CMB data studies. In this 
paper, we concentrate on two particular components: the CMB and the SZ components. 
For such processes, state-of-the-art blind separation methods used on CMB data are: 

- JADE which is a classical Independent Component Analysis technique based on fourth 
order statistics. Its effectiveness at extracting non- Gaussian components such as the 
SZ map was shown in [10]. 

- Spectral Matching ICA (SMICA) (see [5] and [9]) has been devised to accurately 
separate the CMB component. SMICA assumes the case of mixed stationary Gaussian 
components in a noisy environment. It is based on second order statistics. In the Fourier 
representation, colored stationary Gaussian components are discernible based on the 
diversity of their power spectra. SMICA is then well adapted to Gaussian components 
such as CMB. 

Neither of the aforementioned techniques is able to effectively extract both the SZ and 
CMB maps. In this paper we propose a novel sparsity-based component separation tech- 
nique coined Generalized Morphological Component Analysis (GMCA) which turns out 
to be well suited for the recovery of CMB and SZ components. Section 1 describes the 
GMCA model and the algorithm proposed to solve the corresponding optimization prob- 
lem. Numerical experiments are given which illustrate the astounding performances of 
GMCA for CMB and SZ extraction in Section 2. Finally, we show in Section 2.2 that 
GMCA is versatile enough to account for physical priors. 
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1. Generalized Morphological Component Analysis 

1.1. The GMCA model 

In the previous section, we introduced the hnear mixture model in Equation 1. We fur- 
ther assume that all the protagonists of the model in Equation 1 are random components 
(variables or vectors). More particularly, the entries of the noise matrix N are assumed to 
be independently distributed according to a zero mean Gaussian distribution with vari- 
ance o"? depending on the detector. Prom physical considerations, N models instrumental 
noise the level of which varies independently from one detector to another. N is thus a 
random Gaussian variable with zero mean and covariance matrix Fn = diag(cr^, • • • , ct^). 
In practice, as the detectors are assumed to be accurately calibrated, Fn is known with 
high precision. The log-likelihood function is then the following one : 

logP(X|A, S,Fn) = -^||X - AS||2^r. + C (2) 

where C is a constant. The notation lUlnr stands for the Frobenius norm of Y in 

II 11^, In 

the noise covariance metric : ||1^||2 — Trace (Y-^Fn~^Y). From a Bayesian point of 
view, adding physical priors should help the separation task. We first assume no particu- 
lar knowledge about the emission laws of the components modeled by A. For simplicity, 
we consider that each entry of the mixing matrix A is i.i.dj^from a uniform zero mean 
distribution. Note that it would be possible to add some physical constraint on the emis- 
sion laws reflected in A. 

In the general case, source separation is merely a question of diversity and contrast be- 
tween the sources (see [4]). For instance, on the one hand JADE relies on non-Gaussianity 
to distinguish between the sources. On the other, SMICA takes advantage of the di- 
versity of the mixed components' power spectra to achieve the separation task. "Non- 
Gaussianity" and "power spectra diversity" are contrasts between the sources. A com- 
bination of both characteristics, "Non-Gaussianity" and "power spectra diversity" , was 
also proposed to separate CMB from kinetic SZ signal which are otherwise undistin- 
guishable [7] . Recent work has already emphasized on sparsity as a source of diversity to 
improve component separation (see [14] and [2]). In that setting, each source {sj}j=i,.-- ,ri 
is assumed to be sparse in a representation (potentially overcomplete) V. Formally, V 
is a fixed dictionary of signal waveforms written as a T x t matrix. We define the set 
of projection coefficients Oj such that : Vj G {1, ■ ' ' j'^}) = ctj^- Any source sj is 
said to be sparse in V if most of the entries of aj are nearly zero and only a few have 
"significant" amplitudes. When V is overcomplete (T > t), T> is called a dictionary. The 
attractiveness of overcomplete representations in image processing theory resides in the 
potentially very sparse representations that they make possible (see e.g. [6] and references 
therein) . In the field of basic source separation we showed in [2] that morphological diver- 
sity and sparsity are key properties leading to better separation. We noticed that the gist 
of sparsity-based source separation methods leans on the rationale : ^'independent sources 
are distinctly sparse in a dictionary P" . In that study, we considered the simple case of 
morphologically different sources : components were assumed to be sparsely represented 
in different sub-dictionaries. We illustrated that such a sparsity based prior provides a 

^ Independently and identically distributed. 
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very effective way to distinguish between sources. In the present paper, we focus on a 
more general setting : the sources can have similar morphologies (i.e. all the sources are 
sparsely represented over the whole P). When the overcomplete dictionary P is made of 
the union of D orthonormal bases (i.e. V = ■ ■ ■ ,$d]) then each source is modeled 
as the linear combination of D so-called morphological components (see [11] for details 
on Morphological Component Analysis) - each morphological component being sparse in 
a different orthonormal basis - ,^d}' 

D D 

Vj e {1, • • • , n}, Sj = <pjk = ajk'^k (3) 

k=l k=l 

From a statistical viewpoint, we assume that the entries of Ojk = ipjk^^ are i.i.d. from 
a Laplacian probability distribution with scale parameter 1 //x: 

P{<Pjk) oc exp {-n\\(pjk^l\\i) (4) 

where the ^i-norm ||.||i stands for ||x||i = Ylp=i \^[p\ \ ™ which x[p] is the p-th entry of x. 
In practice, the Laplacian prior is well adapted to model leptokurtic sparse signals. We 
classically assume that the morphological components are statistically mutually indepen- 
dent : P(S) = Ylj k -^(Vjk)- Estimating the sources S is then equivalent to estimating 
the set of morphological components {'i^jk}j=i,---,n;k=i,- - .D- In this Bayesian context, we 
propose to estimate those morphological components {fjk} and the mixing matrix A 
from a maximum a posteriori (MAP) leading to the following optimization problem: 

{{^jk},A} = arg max P(X| A, Tn) J] ^(^i't)^(A) (5) 

where we further assumed that the morphological components {^Pjk} are independent of 
A. Owing to Equations 2 and 4, the mixing matrix A and the morphological components 
{<Pjk} are obtained by minimizing the following negative log a posteriori: 

n D 

{{^jk},A.} = arg min ||X - ASH^^r^ + ^^Z] ll^jfe^fe 111 (6) 
^'^^''J' j=ik=i 

where Vj G {1, • • • ,n}, sj = '^k=i ^jk- Equation 6 leads to the GMCA estimates of 
the sources and the mixing matrix in a general sparse component separation context. 
Interestingly, in the case of CMB data, the sources we look for (CMB, galactic dust and 
SZ) are quite sparse in the same unique orthonormal wavelet basis. The dictionary V 
then reduces to a single orthonormal basis In that case, since $ is unitary, Equation 6 
can be rewritten as follows : 

ja, a| = argmin ||X$^ - AaHj^rN + 2Ai||a||i 

= argmin /^(a. A) = argmin„ a/o(A, a) + 2ju/i(a) (7) 

a, A 

where a = S$"^. Note that the estimation is done in the sparse representation $ requiring 
a single transform of the data X^E"^. To remain computationally efficient, GMCA relies on 
practical transforms which generally involve fast implicit operators (typical complexity 
of O (t) or O (t logt)). In [14], the authors also used a unique orthonormal wavelet basis. 
While a gradient descent is used in [14], we use a fast and efficient iterative thresholding 
optimization scheme which we describe in the next section. 
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1.2. Solving the optimization problem 



The maximum a posteriori estimates of the coefficients cx and the mixing matrix in 
Equation 7 lead to a non-convex minimization problem. Note that in Equation 7 the 
functional to be minimized suffers from several invariances : any permutation or rescal- 
ing of the sources and the mixing matrix leaves the product Act unaltered. The scale 
invariance is computationally alleviated by forcing the columns of A to have unit £2 
norm : Vi £ {1, • • • , n}, a* a* = 1 where a* is the i-th column of A. 
As solutions of problem (7) have no explicit formulation, we propose solving it by means 
of a block-coordinate relaxation iterative algorithm such that each iteration (h) is de- 
composed into two steps : (i) estimation of the sources S assuming the mixing matrix is 
fixed to its current estimate A^''"^) and (ii) estimation of the mixing matrix assuming 
the sources are fixed to their current estimates S^''^ It is not difficult to see that the 
objective MAP functional in (7) is continuous on its effective domain and has compact 
level sets. Moreover, this objective function is convex in the source coefficient vectors 
(ai, . . . , an), and /o has an open domain, is continuous and Gateaux diiferentiable. Thus 
by [13, Theorem 4.1], the iterates generated by our alternating algorithm are defined and 
bounded, and each accumulation point is a stationary point of the MAP functional. In 
other words, our iterative algorithm will converge. Hence, at iteration (h), the sources 
are estimated from a maximum a posteriori assuming A = A^''"^^ By classical ideas 
in convex analysis, a necessary condition for a to be a minimizer is that the zero is an 
element of the subdifferential of the objective at ct. We calculat^^ 

A) = -2A^rN"^(X$^ - Aa) + 2/i5c||a||i (8) 

where 9ct||cK||i is defined as (owing to the separability of the prior): 

Uj,k = sign(aj_fc), 0^,^/0 



dMU = {Ue 



hUXt 



Uj,k e[-l,l], aj,fe = 



Hence, Equation 8 can be rewritten equivalently as two conditions leading to the following 
(proximal) fixed point equation: 



aj,k 



0, if |(A^rN-'XcD^)._J <^ 



A^Pn ^(X^>^ - Aa) = /i sign (a) , otherwise. 



(9) 



Unfortunately, Equation 9 has no closed-form solution in general. It must be iterated and 
is thus computationally demanding. Fortunately, it can be simplified when A has nearly 
orthogonal columns in the noise covariance matrix {i.e. A-^Pn~^ A ~ diag ( A-^Pn~^ A ) ). 



Let C = ("a^Pn ^a) A^Pn ^X$^, Equation 9 boils down to the following set of 



equations Vj £ {!,••• , n}: 



d,.fc =0, if |C,- ,1 < M^'^)a| 

(10) 



ctj = [C]j — fiaj sign(dj) , otherwise. 
For clarity, we drop the superscript (h — 1) and write A = AC'"^). 
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where [Y]j is the j-th row of Y. In practice, even if the approximation we make is not 
strictly vahd, such a simpHfication leads to good computational results. These equa- 
tions are known as soft-thresholding with threshold We define ST5(.), the soft- 
thresholding operator with threshold d. At iteration (h), the sources are thus estimated 
such that: 

df ) = ST^(.).| ([C] .) (11) 

The jth source is reconstructed as Sj^^ = a^^^^. The mixing matrix A is then estimated 
by a maximum likelihood estimate amounting to a simple least-squares update assuming 
S is fixed. The GMCA algorithm is then described as follows : 



1. Set the number of iterations /max and thresholds 5^°^ = 

2. While each /iC') is higher than a given lower bound fJ-min (s g- can depend on the noise variance), 
— Proceed with the following iteration to estimate source coefficients cx at iteration h assuming A is 



fixed: a('^) = ST^(.)^ 



Update A assuming ct is fixed : A*^'') = X$^a^ (aa^) 
- Decrease the threshold /i^'') following a given strategy 



Note that the overall optimization scheme is based on an iterative and alternate thresh- 
olding algorithm involving a coarse to fine estimation process. Indeed, coarse versions of 
the sources (i.e. containing the most "significant" features of the sources) are first com- 
puted with high values of fi'^^K In the early stages of the algorithm, the mixing matrix 
is then estimated from the most "significant" features of the sources which are less per- 
turbed by noise. The estimation of A and S is then refined at each iteration as fi^'^^ (and 
thus the thresholds {M*^''^<^|}i=i,--- ,«) decreases towards a final value i^min- We already 
used this minimization scheme in [2] where this optimization process provided robustness 
and helped convergence even in a noisy context. Experiments in Section 2 illustrate that 
it achieves good results with GMCA as well. 



2. Application to CMB and SZ reconstruction 

2.1. Blind component separation 

The method described above was applied to synthetic data composed of m = 6 mix- 
tures of n = 3 sources : CMB, galactic dust emission and SZ maps illustrated in Fig- 
ure 1 and 2. Following [5] and [9], we do not take into account the emission at high 
frequency from the fluctuations of infra-red galaxies. The synthetic data mimic the ob- 
servations that will be acquired in the six frequency channels of Planck-HFI namely : 
100,143,217,353,545 and 857 GHz, as shown on Figure 2. White Gaussian noise N is 
added with diagonal covariance matrix Fn reflecting the foreseen Planck-HFI noise levels. 
Experiments were led with 7 global noise levels with SNR from 1.7 to 16.7dB such that the 
experimental noise covariance Fn was proportional to the nominal noise covariance. Note 
that the nominal Planck-HFI global noise level is about lOdB. Each measurement point 
was computed from 30 experiments involving random noise, randomly chosen sources 
from a data set of several simulated CMB, galactic dust and SZ 256 x 256 maps. The 
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Fig. 1. The simulated sources - Left: CMB. Middle: galactic dust emission. Right: SZ map. 




Fig. 2. The observed CMB data - global SNR = 2.7dB 



astrophysical components and the mixture maps were generated as in [9] according to 
equation (1) based on model or experimental emission laws, possibly extrapolated, of the 
individual components. Separation was obtained with GMCA using a single orthonor- 
mal wavelet basis. Figure 3 depicts the average correlation coefficients over experiments 
between the estimated source maps and the true source maps. Figure 3 upper left panel 
shows the correlation coefficient between the true simulated CMB map and the one es- 
timated by JADE {dotted line with □), SMICA {dashed line with o) and GMCA {solid 
line). The CMB map is well estimated by SMICA, which indeed was designed for the 
blind separation of stationary colored Gaussian processes, but not as well using JADE as 
one might have expected. GMCA turns out to perform similarly to SMICA. In the second 
line on the left of Figure 3, galactic dust is well estimated by both GMCA and SMICA. 
The SMICA estimates seem to have a slightly higher variance than GMCA estimates for 
higher global noise levels (SNR lower than 5 dB). Finally, the picture in the third line on 
the left shows that GMCA gives better estimates of the SZ map than SMICA when the 
noise variance increases. The right panels provide the dispersion {i.e. standard deviation) 
of the correlation coefficients of the sources estimates. It appears that GMCA is a gen- 
eral method yielding simultaneous SZ and CMB estimates comparable to state-of-the-art 
blind separation techniques which seem mostly dedicated to individual components. 
In a noisy context, assessing separation techniques turns out to be more accurate using 
a mixing matrix criterion. We define the mixing matrix criterion Aa = ||I — PA~^A||i^i 
(where P is a matrix that reduces the scale/permutation indeterminacy of the mixing 
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model, and ||.||i^i is the entry wise £i matrix norm). Indeed, when A is perfectly estimated, 
it is equal to A up to scale and permutation. As we entirely manage our experiments, 
the true sources and mixing matrix are known and thus P can be computed easily. The 
mixing matrix criterion is thus strictly positive unless the mixing matrix is perfectly 
estimated up to scale and permutation. This mixing matrix criterion is experimentally 
much more sensitive to separation error. The bottom right panel of Figure 3 illustrates 
the behavior of the mixing matrix criterion Aa with JADE, SMICA and GMCA as the 
global noise variance varies. GMCA clearly outperforms SMICA and JADE when applied 
to CMB data. 

2.2. Adding some physical constraint : the versatility of GMCA 

In practice, the separation task is only partly blind. Indeed, the CMB emission law is 
extremely well-known. In this section, we illustrate that GMCA is versatile enough to 
account for such prior knowledge. In the following experiment, CMB-GMCA has been 
designed by constraining the column of the mixing matrix A related to CMB to its 
true value. This is equivalent to placing a strict prior on the CMB column of A; that is 
p(^(icmb^ _ g(^Qcm.b _ ^cmb^ where 6{.) is the Dirac measure and Og™''' is the true simulated 
CMB emission law in the frequency range of Planck-HFI. Figure 4 shows the correlation 
coelficients between the true source maps and the source maps estimated using GMCA 
with and without the CMB prior. As expected, the top left picture of Figure 4 shows 
that assuming a^^^ is known improves the estimation of CMB. Interestingly, the galactic 
dust map (in the middle on the left of Figure 4) is also better estimated. Furthermore, 
the CMB-GMCA SZ map estimate is likely to have a lower variance (bottom-right panel 
of Figure 4) . Moreover, it is likely to provide more robustness to the SZ and galactic dust 
estimates thus enhancing the global separation performances. 

3. Conclusion 

In this paper we underlined that recovering information from CMB data requires solv- 
ing a blind source separation issue (BSS). Several BSS techniques have already been 
applied to CMB data without providing good global performances i.e. on all components 
simultaneously. In this paper, we provide a sparsity-based source separation method 
coined Generalized Morphological Component Analysis (GMCA) which turns to give 
astounding results to effectively recover both CMB and SZ maps. In that context, spar- 
sity enhances the contrast between the sources leading to an improved separation task 
even in a noisy context. In the blind case, when no prior knowledge is assumed on the 
emission laws of the components, GMCA outperforms state-of-the-art blind component 
separation techniques already applied to CMB data. Furthermore, GMCA is versatile 
enough to easily include some prior knowledge of the emission laws of the components. 
This is an extremely valuable feature of the proposed method in the case of CMB data 
analysis. Indeed, the CMB has a known black-body spectrum. Including this information 
in the GMCA algorithm enhances the source separation globally. In the present work, we 
have chosen not to include prior knowledge on the SZ spectral signature. Adding such 
prior would lead to even better separation results. Future work will be devoted to taking 
advantage of GMCA's versatility to adapt to more complex physical models. 
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Fig. 3. Left Column : Mean value of the correlation coefficients between the estimated 
source map and the true source map - Right Column : Dispersion of these correlation 
coefficients : First line : CMB. Second line: galactic dust. Third line: SZ map. Fourth line: 

mixing matrix criterion Aa Legend : JADE : dotted line with □ - SMICA : dashed line with o - GMCA 
: solid line. Abscissa : SNR in dB. 
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Fig. 4. Left Column : Mean value of the correlation coefficients between the estimated 
source map and the true source map - Right Column : Dispersion of these correlation 
coefficients : First line : CMB. Second line: galactic dust. Third line: SZ map. Legend : GMCA 
assuming that the CMB emission law is known : dotted line - GMCA : solid line. Abscissa : SNR in 
dB. 
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